# ----------------------------------------------------------------------------------------------------------------------
# Replication script for Tables 2and 3 in van Ditmars and Ksiazkiewicz
# "The gender gap in political interest: heritability, gendered political socialization, and the enriched environment hypothesis"
# Politics and the Life Sciences

# This script was adapted by Aleks Ksiazkiewicz from the scripts below from the International Statistical Genetics Workshop at UC Boulder in 2018.

#   Based on: oneACE5c.R
#   Author: Eveline de Zeeuw
#   Date: 04 03 2018
#   which was adapted from: oneADE5ca.R
#   Author: Hermine Maes
#	  Date: 01 04 2018

# Note that this script requires the version of OpenMx available here: source('https://vipbg.vcu.edu/vipbg/OpenMx2/software/getOpenMx.R')
# It will not run on the version of OpenMx available on CRAN because it requires NPSOL
# Details here: https://openmx.ssri.psu.edu/installing-openmx

# This script also requires a helper function script to display the results of the analysis (miFunctionsageandsex.R)

# Place both of these scripts in the folder with the data to be analyzed (twinlife.csv) and set the working directory accordingly.

# The data for this script are v 3-0-0 of the TwinLife dataset cleaned with the cleaning script provided for Tables 2 and 3.

# ----------------------------------------------------------------------------------------------------------------------

setwd("")

# Remove all Objects
rm(list=ls(all=TRUE))

# Load Libraries & Options 
library(OpenMx)
library(psych)
source("miFunctionsageandsex.R")
mxOption( NULL, "Default optimizer", "NPSOL" )

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE DATA

# Load Data
Data <- read.table ('twinlife.csv', header=T, sep=',', na.strings=".")
describe(Data, skew = F)
dim(Data)
head(Data)

varnames <- rbind("interest_")

# Select Variables for Analysis
vars      <- c(varnames[1,1])          # list of variables names
nv        <- 1                         # number of variables
ntv       <- nv*2                      # number of total variables
selVars   <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")

# Select Data for Analysis
mz1Data	<- subset(Data, (mzage==1 & mzsex==1), selVars) # 11-12 yo male mz
dz1Data	<- subset(Data, (mzage==2 & mzsex==2), selVars) # 11-12 yo male dz
mz2Data	<- subset(Data, (mzage==3 & mzsex==1), selVars) # 17-18 yo male mz
dz2Data	<- subset(Data, (mzage==4 & mzsex==2), selVars) # 17-18 yo male dz
mz3Data	<- subset(Data, (mzage==5 & mzsex==1), selVars) # 22-25 yo male mz
dz3Data	<- subset(Data, (mzage==6 & mzsex==2), selVars) # 22-25 yo male dz
mz4Data	<- subset(Data, (mzage==1 & mzsex==3), selVars) # 11-12 yo female mz
dz4Data	<- subset(Data, (mzage==2 & mzsex==4), selVars) # 11-12 yo female dz
mz5Data	<- subset(Data, (mzage==3 & mzsex==3), selVars) # 17-18 yo female mz
dz5Data	<- subset(Data, (mzage==4 & mzsex==4), selVars) # 17-18 yo female dz
mz6Data	<- subset(Data, (mzage==5 & mzsex==3), selVars) # 22-25 yo female mz
dz6Data	<- subset(Data, (mzage==6 & mzsex==4), selVars) # 22-25 yo female dz



# Generate Descriptive Statistics 
colMeans(mz1Data,na.rm=TRUE) 
colMeans(dz1Data,na.rm=TRUE) 
colMeans(mz2Data,na.rm=TRUE) 
colMeans(dz2Data,na.rm=TRUE) 
colMeans(mz3Data,na.rm=TRUE) 
colMeans(dz3Data,na.rm=TRUE) 
colMeans(mz4Data,na.rm=TRUE) 
colMeans(dz4Data,na.rm=TRUE) 
colMeans(mz5Data,na.rm=TRUE) 
colMeans(dz5Data,na.rm=TRUE) 
colMeans(mz6Data,na.rm=TRUE) 
colMeans(dz6Data,na.rm=TRUE) 

cov(mz1Data,use="complete") 
cov(dz1Data,use="complete")
cov(mz2Data,use="complete") 
cov(dz2Data,use="complete")
cov(mz3Data,use="complete") 
cov(dz3Data,use="complete")
cov(mz4Data,use="complete") 
cov(dz4Data,use="complete")
cov(mz5Data,use="complete") 
cov(dz5Data,use="complete")
cov(mz6Data,use="complete") 
cov(dz6Data,use="complete")

cor(mz1Data,use="complete") 
cor(dz1Data,use="complete")
cor(mz2Data,use="complete") 
cor(dz2Data,use="complete")
cor(mz3Data,use="complete") 
cor(dz3Data,use="complete")
cor(mz4Data,use="complete") 
cor(dz4Data,use="complete")
cor(mz5Data,use="complete") 
cor(dz5Data,use="complete")
cor(mz6Data,use="complete") 
cor(dz6Data,use="complete")

# ----------------------------------------------------------------------------------------------------------------------
# PREPARE MODEL

# Set Starting Values
svMe1	<- .05	  # start value for means for age1
svMe2	<- -.05  	# start value for means for age2
svMe3	<- -.05  	# start value for means for age3
svMe4	<- .05	  # start value for means for age1
svMe5 <- -.05  	# start value for means for age2
svMe6	<- -.05  	# start value for means for age3
svPa1	<- .8	# start value for a for age1
svPc1	<- .2	# start value for c for age1
svPe1	<- .3	# start value for e for age1
svPa2	<- .8	# start value for a for age2 
svPc2	<- .2	# start value for c for age2
svPe2	<- .3	# start value for e for age2
svPa3	<- .8	# start value for a for age3 
svPc3	<- .2	# start value for c for age3
svPe3	<- .3	# start value for e for age3
svPa4	<- .8	# start value for a for age1
svPc4	<- .2	# start value for c for age1
svPe4	<- .3	# start value for e for age1
svPa5	<- .8	# start value for a for age2 
svPc5	<- .2	# start value for c for age2
svPe5	<- .3	# start value for e for age2
svPa6	<- .8	# start value for a for age3 
svPc6	<- .2	# start value for c for age3
svPe6	<- .3	# start value for e for age3


# Create Algebra for expected Mean Matrices
meanG1     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svMe1, labels="mean1", name="meanG1" )
meanG2     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svMe2, labels="mean2", name="meanG2" )
meanG3     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svMe3, labels="mean3", name="meanG3" )
meanG4     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svMe4, labels="mean4", name="meanG4" )
meanG5     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svMe5, labels="mean5", name="meanG5" )
meanG6     <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, values=svMe6, labels="mean6", name="meanG6" )
expMeanMZ1    <- mxAlgebra( expression= cbind(meanG1, meanG1), name="expMeanMZ1" )
expMeanDZ1    <- mxAlgebra( expression= cbind(meanG1, meanG1), name="expMeanDZ1" )
expMeanMZ2    <- mxAlgebra( expression= cbind(meanG2, meanG2), name="expMeanMZ2" )
expMeanDZ2    <- mxAlgebra( expression= cbind(meanG2, meanG2), name="expMeanDZ2" )
expMeanMZ3    <- mxAlgebra( expression= cbind(meanG3, meanG3), name="expMeanMZ3" )
expMeanDZ3    <- mxAlgebra( expression= cbind(meanG3, meanG3), name="expMeanDZ3" )
expMeanMZ4    <- mxAlgebra( expression= cbind(meanG4, meanG4), name="expMeanMZ4" )
expMeanDZ4    <- mxAlgebra( expression= cbind(meanG4, meanG4), name="expMeanDZ4" )
expMeanMZ5    <- mxAlgebra( expression= cbind(meanG5, meanG5), name="expMeanMZ5" )
expMeanDZ5    <- mxAlgebra( expression= cbind(meanG5, meanG5), name="expMeanDZ5" )
expMeanMZ6    <- mxAlgebra( expression= cbind(meanG6, meanG6), name="expMeanMZ6" )
expMeanDZ6    <- mxAlgebra( expression= cbind(meanG6, meanG6), name="expMeanDZ6" )

# Create Matrices for Path Coefficients
pathA1     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa1, labels="a111", name="a1" ) 
pathC1     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPc1, labels="c111", name="c1" ) 
pathE1     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe1, labels="e111", name="e1" )
pathA2     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa2, labels="a211", name="a2" ) 
pathC2     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPc2, labels="c211", name="c2" ) 
pathE2     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe2, labels="e211", name="e2" )
pathA3     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa3, labels="a311", name="a3" ) 
pathC3     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPc3, labels="c311", name="c3" ) 
pathE3     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe3, labels="e311", name="e3" )
pathA4     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa4, labels="a411", name="a4" ) 
pathC4     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPc4, labels="c411", name="c4" ) 
pathE4     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe4, labels="e411", name="e4" )
pathA5     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa5, labels="a511", name="a5" ) 
pathC5     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPc5, labels="c511", name="c5" ) 
pathE5     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe5, labels="e511", name="e5" )
pathA6     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPa6, labels="a611", name="a6" ) 
pathC6     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPc6, labels="c611", name="c6" ) 
pathE6     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPe6, labels="e611", name="e6" )

# Create Algebra for Variance Components
covA1      <- mxAlgebra( expression=a1 %*% t(a1), name="A1" ) 
covC1      <- mxAlgebra( expression=c1 %*% t(c1), name="C1" ) 
covE1      <- mxAlgebra( expression=e1 %*% t(e1), name="E1" )
covA2      <- mxAlgebra( expression=a2 %*% t(a2), name="A2" ) 
covC2      <- mxAlgebra( expression=c2 %*% t(c2), name="C2" ) 
covE2      <- mxAlgebra( expression=e2 %*% t(e2), name="E2" )
covA3      <- mxAlgebra( expression=a3 %*% t(a3), name="A3" ) 
covC3      <- mxAlgebra( expression=c3 %*% t(c3), name="C3" ) 
covE3      <- mxAlgebra( expression=e3 %*% t(e3), name="E3" )
covA4      <- mxAlgebra( expression=a4 %*% t(a4), name="A4" ) 
covC4      <- mxAlgebra( expression=c4 %*% t(c4), name="C4" ) 
covE4      <- mxAlgebra( expression=e4 %*% t(e4), name="E4" )
covA5      <- mxAlgebra( expression=a5 %*% t(a5), name="A5" ) 
covC5      <- mxAlgebra( expression=c5 %*% t(c5), name="C5" ) 
covE5      <- mxAlgebra( expression=e5 %*% t(e5), name="E5" )
covA6      <- mxAlgebra( expression=a6 %*% t(a6), name="A6" ) 
covC6      <- mxAlgebra( expression=c6 %*% t(c6), name="C6" ) 
covE6      <- mxAlgebra( expression=e6 %*% t(e6), name="E6" )

# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP1     <- mxAlgebra( expression= A1+C1+E1, name="V1" )
covP2     <- mxAlgebra( expression= A2+C2+E2, name="V2" )
covP3     <- mxAlgebra( expression= A3+C3+E3, name="V3" )
covP4     <- mxAlgebra( expression= A4+C4+E4, name="V4" )
covP5     <- mxAlgebra( expression= A5+C5+E5, name="V5" )
covP6     <- mxAlgebra( expression= A6+C6+E6, name="V6" )

covMZ1    <- mxAlgebra( expression= A1+C1, name="cMZ1" )
covDZ1    <- mxAlgebra( expression= 0.5%x%A1+ C1, name="cDZ1" )
covMZ2    <- mxAlgebra( expression= A2+C2, name="cMZ2" )
covDZ2    <- mxAlgebra( expression= 0.5%x%A2+ C2, name="cDZ2" )
covMZ3    <- mxAlgebra( expression= A3+C3, name="cMZ3" )
covDZ3    <- mxAlgebra( expression= 0.5%x%A3+ C3, name="cDZ3" )
covMZ4    <- mxAlgebra( expression= A4+C4, name="cMZ4" )
covDZ4    <- mxAlgebra( expression= 0.5%x%A4+ C4, name="cDZ4" )
covMZ5    <- mxAlgebra( expression= A5+C5, name="cMZ5" )
covDZ5    <- mxAlgebra( expression= 0.5%x%A5+ C5, name="cDZ5" )
covMZ6    <- mxAlgebra( expression= A6+C6, name="cMZ6" )
covDZ6    <- mxAlgebra( expression= 0.5%x%A6+ C6, name="cDZ6" )

expCovMZ1 <- mxAlgebra( expression= rbind( cbind(V1, cMZ1), cbind(t(cMZ1), V1)), name="expCovMZ1" )
expCovDZ1 <- mxAlgebra( expression= rbind( cbind(V1, cDZ1), cbind(t(cDZ1), V1)), name="expCovDZ1" )
expCovMZ2 <- mxAlgebra( expression= rbind( cbind(V2, cMZ2), cbind(t(cMZ2), V2)), name="expCovMZ2" )
expCovDZ2 <- mxAlgebra( expression= rbind( cbind(V2, cDZ2), cbind(t(cDZ2), V2)), name="expCovDZ2" )
expCovMZ3 <- mxAlgebra( expression= rbind( cbind(V3, cMZ3), cbind(t(cMZ3), V3)), name="expCovMZ3" )
expCovDZ3 <- mxAlgebra( expression= rbind( cbind(V3, cDZ3), cbind(t(cDZ3), V3)), name="expCovDZ3" )
expCovMZ4 <- mxAlgebra( expression= rbind( cbind(V4, cMZ4), cbind(t(cMZ4), V4)), name="expCovMZ4" )
expCovDZ4 <- mxAlgebra( expression= rbind( cbind(V4, cDZ4), cbind(t(cDZ4), V4)), name="expCovDZ4" )
expCovMZ5 <- mxAlgebra( expression= rbind( cbind(V5, cMZ5), cbind(t(cMZ5), V5)), name="expCovMZ5" )
expCovDZ5 <- mxAlgebra( expression= rbind( cbind(V5, cDZ5), cbind(t(cDZ5), V5)), name="expCovDZ5" )
expCovMZ6 <- mxAlgebra( expression= rbind( cbind(V6, cMZ6), cbind(t(cMZ6), V6)), name="expCovMZ6" )
expCovDZ6 <- mxAlgebra( expression= rbind( cbind(V6, cDZ6), cbind(t(cDZ6), V6)), name="expCovDZ6" )


# Create Data Objects for Multiple Groups
dataMZ1    <- mxData( observed=mz1Data, type="raw" )
dataDZ1    <- mxData( observed=dz1Data, type="raw" )
dataMZ2    <- mxData( observed=mz2Data, type="raw" )
dataDZ2    <- mxData( observed=dz2Data, type="raw" )
dataMZ3    <- mxData( observed=mz3Data, type="raw" )
dataDZ3    <- mxData( observed=dz3Data, type="raw" )
dataMZ4    <- mxData( observed=mz4Data, type="raw" )
dataDZ4    <- mxData( observed=dz4Data, type="raw" )
dataMZ5    <- mxData( observed=mz5Data, type="raw" )
dataDZ5    <- mxData( observed=dz5Data, type="raw" )
dataMZ6    <- mxData( observed=mz6Data, type="raw" )
dataDZ6    <- mxData( observed=dz6Data, type="raw" )

# Create Expectation Objects for Multiple Groups
expMZ1     <- mxExpectationNormal( covariance="expCovMZ1", means="expMeanMZ1", dimnames=selVars ) 
expDZ1     <- mxExpectationNormal( covariance="expCovDZ1", means="expMeanDZ1", dimnames=selVars ) 
expMZ2     <- mxExpectationNormal( covariance="expCovMZ2", means="expMeanMZ2", dimnames=selVars ) 
expDZ2     <- mxExpectationNormal( covariance="expCovDZ2", means="expMeanDZ2", dimnames=selVars ) 
expMZ3     <- mxExpectationNormal( covariance="expCovMZ3", means="expMeanMZ3", dimnames=selVars ) 
expDZ3     <- mxExpectationNormal( covariance="expCovDZ3", means="expMeanDZ3", dimnames=selVars ) 
expMZ4     <- mxExpectationNormal( covariance="expCovMZ4", means="expMeanMZ4", dimnames=selVars ) 
expDZ4     <- mxExpectationNormal( covariance="expCovDZ4", means="expMeanDZ4", dimnames=selVars )
expMZ5     <- mxExpectationNormal( covariance="expCovMZ5", means="expMeanMZ5", dimnames=selVars ) 
expDZ5     <- mxExpectationNormal( covariance="expCovDZ5", means="expMeanDZ5", dimnames=selVars )
expMZ6     <- mxExpectationNormal( covariance="expCovMZ6", means="expMeanMZ6", dimnames=selVars ) 
expDZ6     <- mxExpectationNormal( covariance="expCovDZ6", means="expMeanDZ6", dimnames=selVars )
funML      <- mxFitFunctionML()

# Create Model Objects for Multiple Groups
pars1      <- list( meanG1, pathA1, pathC1, pathE1, covA1, covC1, covE1, covP1 )
pars2      <- list( meanG2, pathA2, pathC2, pathE2, covA2, covC2, covE2, covP2 )
pars3      <- list( meanG3, pathA3, pathC3, pathE3, covA3, covC3, covE3, covP3 )
pars4      <- list( meanG4, pathA4, pathC4, pathE4, covA4, covC4, covE4, covP4 )
pars5      <- list( meanG5, pathA5, pathC5, pathE5, covA5, covC5, covE5, covP5 )
pars6      <- list( meanG6, pathA6, pathC6, pathE6, covA6, covC6, covE6, covP6 )

modelMZ1   <- mxModel( pars1, covMZ1, expCovMZ1, expMeanMZ1, dataMZ1, expMZ1, funML, name="MZ1" ) 
modelDZ1   <- mxModel( pars1, covDZ1, expCovDZ1, expMeanDZ1,  dataDZ1, expDZ1, funML, name="DZ1" ) 
modelMZ2   <- mxModel( pars2, covMZ2, expCovMZ2, expMeanMZ2,  dataMZ2, expMZ2, funML, name="MZ2" ) 
modelDZ2   <- mxModel( pars2, covDZ2, expCovDZ2, expMeanDZ2,  dataDZ2, expDZ2, funML, name="DZ2" ) 
modelMZ3   <- mxModel( pars3, covMZ3, expCovMZ3, expMeanMZ3,  dataMZ3, expMZ3, funML, name="MZ3" ) 
modelDZ3   <- mxModel( pars3, covDZ3, expCovDZ3, expMeanDZ3,  dataDZ3, expDZ3, funML, name="DZ3" ) 
modelMZ4   <- mxModel( pars4, covMZ4, expCovMZ4, expMeanMZ4,  dataMZ4, expMZ4, funML, name="MZ4" ) 
modelDZ4   <- mxModel( pars4, covDZ4, expCovDZ4, expMeanDZ4,  dataDZ4, expDZ4, funML, name="DZ4" ) 
modelMZ5   <- mxModel( pars5, covMZ5, expCovMZ5, expMeanMZ5,  dataMZ5, expMZ5, funML, name="MZ5" ) 
modelDZ5   <- mxModel( pars5, covDZ5, expCovDZ5, expMeanDZ5,  dataDZ5, expDZ5, funML, name="DZ5" ) 
modelMZ6   <- mxModel( pars6, covMZ6, expCovMZ6, expMeanMZ6,  dataMZ6, expMZ6, funML, name="MZ6" ) 
modelDZ6   <- mxModel( pars6, covDZ6, expCovDZ6, expMeanDZ6,  dataDZ6, expDZ6, funML, name="DZ6" ) 

multi     <- mxFitFunctionMultigroup( c("MZ1","DZ1","MZ2","DZ2","MZ3","DZ3","MZ4","DZ4","MZ5","DZ5","MZ6","DZ6") )

# Create Algebra for Variance Components
rowVC     <- rep('VC',nv)
colVC0    <- rep(c('V1','V2','V3','V4','V5','V6'),each=nv)
colVC1    <- rep(c('A1','C1','E1','SA1','SC1','SE1'),each=nv)
colVC2    <- rep(c('A2','C2','E2','SA2','SC2','SE2'),each=nv)
colVC3    <- rep(c('A3','C3','E3','SA3','SC3','SE3'),each=nv)
colVC4    <- rep(c('A4','C4','E4','SA4','SC4','SE4'),each=nv)
colVC5    <- rep(c('A5','C5','E5','SA5','SC5','SE5'),each=nv)
colVC6    <- rep(c('A6','C6','E6','SA6','SC6','SE6'),each=nv)

estVC0    <- mxAlgebra( expression=cbind(V1,V2,V3,V4,V5,V6), name="VC0", dimnames=list(rowVC,colVC0))
estVC1    <- mxAlgebra( expression=cbind(A1,C1,E1,A1/V1,C1/V1,E1/V1), name="VC1", dimnames=list(rowVC,colVC1))
estVC2    <- mxAlgebra( expression=cbind(A2,C2,E2,A2/V2,C2/V2,E2/V2), name="VC2", dimnames=list(rowVC,colVC2))
estVC3    <- mxAlgebra( expression=cbind(A3,C3,E3,A3/V3,C3/V3,E3/V3), name="VC3", dimnames=list(rowVC,colVC3))
estVC4    <- mxAlgebra( expression=cbind(A4,C4,E4,A4/V4,C4/V4,E4/V4), name="VC4", dimnames=list(rowVC,colVC4))
estVC5    <- mxAlgebra( expression=cbind(A5,C5,E5,A5/V5,C5/V5,E5/V5), name="VC5", dimnames=list(rowVC,colVC5))
estVC6    <- mxAlgebra( expression=cbind(A6,C6,E6,A6/V6,C6/V6,E6/V6), name="VC6", dimnames=list(rowVC,colVC6))

# Confidence intervals for A component
rowdiff <- rep('Adiff',nv)
colAdiff <- rep(c('ad12','ad13','ad14','ad23','ad25','ad36','ad45','ad46','ad56'),each=nv)
estAdiff <- mxAlgebra( expression=cbind(A1/V1-A2/V2,A1/V1-A3/V3,A1/V1-A4/V4,A2/V2-A3/V3,A2/V2-A5/V5,A3/V3-A6/V6,A4/V4-A5/V5,A4/V4-A6/V6,A5/V5-A6/V6), name="Adiff",dimnames=list(rowdiff,colAdiff))
estAdiff90 <- mxAlgebra( expression=cbind(A1/V1-A2/V2,A1/V1-A3/V3,A1/V1-A4/V4,A2/V2-A3/V3,A2/V2-A5/V5,A3/V3-A6/V6,A4/V4-A5/V5,A4/V4-A6/V6,A5/V5-A6/V6), name="Adiff90",dimnames=list(rowdiff,colAdiff))
colCdiff <- rep(c('cd12','cd13','cd14','cd23','cd25','cd36','cd45','cd46','cd56'),each=nv)
estCdiff <- mxAlgebra( expression=cbind(C1/V1-C2/V2,C1/V1-C3/V3,C1/V1-C4/V4,C2/V2-C3/V3,C2/V2-C5/V5,C3/V3-C6/V6,C4/V4-C5/V5,C4/V4-C6/V6,C5/V5-C6/V6), name="Cdiff",dimnames=list(rowdiff,colCdiff))
estCdiff90 <- mxAlgebra( expression=cbind(C1/V1-C2/V2,C1/V1-C3/V3,C1/V1-C4/V4,C2/V2-C3/V3,C2/V2-C5/V5,C3/V3-C6/V6,C4/V4-C5/V5,C4/V4-C6/V6,C5/V5-C6/V6), name="Cdiff90",dimnames=list(rowdiff,colCdiff))
colEdiff <- rep(c('ed12','ed13','ed14','ed23','ed25','ed36','ed45','ed46','ed56'),each=nv)
estEdiff <- mxAlgebra( expression=cbind(E1/V1-E2/V2,E1/V1-E3/V3,E1/V1-E4/V4,E2/V2-E3/V3,E2/V2-E5/V5,E3/V3-E6/V6,E4/V4-E5/V5,E4/V4-E6/V6,E5/V5-E6/V6), name="Ediff",dimnames=list(rowdiff,colEdiff))
estEdiff90 <- mxAlgebra( expression=cbind(E1/V1-E2/V2,E1/V1-E3/V3,E1/V1-E4/V4,E2/V2-E3/V3,E2/V2-E5/V5,E3/V3-E6/V6,E4/V4-E5/V5,E4/V4-E6/V6,E5/V5-E6/V6), name="Ediff90",dimnames=list(rowdiff,colEdiff))
colVdiff <- rep(c('vd12','vd13','vd14','vd23','vd25','vd36','vd45','vd46','vd56'),each=nv)
estVdiff <- mxAlgebra( expression=cbind(V1-V2,V1-V3,V1-V4,V2-V3,V2-V5,V3-V6,V4-V5,V4-V6,V5-V6), name="Vdiff",dimnames=list(rowdiff,colAdiff))
estVdiff90 <- mxAlgebra( expression=cbind(V1-V2,V1-V3,V1-V4,V2-V3,V2-V5,V3-V6,V4-V5,V4-V6,V5-V6), name="Vdiff90",dimnames=list(rowdiff,colAdiff))

# Create Confidence Interval Objects
ciACE     <- mxCI( c("VC1[1,4:6]","VC2[1,4:6]","VC3[1,4:6]","VC4[1,4:6]","VC5[1,4:6]","VC6[1,4:6]","VC0[1,1:6]") )
ciAdiff     <- mxCI( c("Adiff[1,1:9]") )
ciAdiff90     <- mxCI( c("Adiff90[1,1:9]"), interval = 0.90 )
ciCdiff     <- mxCI( c("Cdiff[1,1:9]") )
ciCdiff90     <- mxCI( c("Cdiff90[1,1:9]"), interval = 0.90 )
ciEdiff     <- mxCI( c("Ediff[1,1:9]") )
ciEdiff90     <- mxCI( c("Ediff90[1,1:9]"), interval = 0.90 )
ciVdiff     <- mxCI( c("Vdiff[1,1:9]") )
ciVdiff90     <- mxCI( c("Vdiff90[1,1:9]"), interval = 0.90 )

# Build Model with Confidence Intervals
modelACErq  <- mxModel( "oneACErq5c", pars1, pars2, pars3, pars4, pars5, pars6, 
                        modelMZ1, modelDZ1, modelMZ2, modelDZ2, modelMZ3, modelDZ3, modelMZ4, modelDZ4, modelMZ5, modelDZ5, modelMZ6, modelDZ6, 
                        multi, 
                        estVC1, estVC2, estVC3, estVC4, estVC5, estVC6, estAdiff, estAdiff90, estCdiff, estCdiff90, estEdiff, estEdiff90, estVC0, estVdiff, estVdiff90,
                        ciACE, ciAdiff, ciAdiff90, ciCdiff, ciCdiff90, ciEdiff, ciEdiff90, ciVdiff, ciVdiff90)

# ----------------------------------------------------------------------------------------------------------------------
# RUN SUBMODELS

# Run Scalar Sex-Limitation ACE model
modelACEq <- mxModel( modelACErq, name="oneACEq5c" )
# modelACEq <- omxSetParameters( modelACEq, labels="ra11", free=FALSE, values=1 )
fitACEq   <- mxRun( modelACEq, intervals=T )
sumACEq   <- summary( fitACEq )
# mxCompare( fitACErq, fitACEq )
fitGofs(fitACEq)
fitEstVCageandsex(fitACEq)



# Table 2 - column 1 (All)
# Run No age or sex differences model
modelACEnone <- mxModel( modelACEq, name="oneACE5none" )
modelACEnone <- omxSetParameters( modelACEnone, labels="a111", newlabels="a11", values=.6 )
modelACEnone <- omxSetParameters( modelACEnone, labels="a211", newlabels="a11", values=.6 )
modelACEnone <- omxSetParameters( modelACEnone, labels="a311", newlabels="a11", values=.6 )
modelACEnone <- omxSetParameters( modelACEnone, labels="a411", newlabels="a11", values=.6 )
modelACEnone <- omxSetParameters( modelACEnone, labels="a511", newlabels="a11", values=.6 )
modelACEnone <- omxSetParameters( modelACEnone, labels="a611", newlabels="a11", values=.6 )
modelACEnone <- omxSetParameters( modelACEnone, labels="c111", newlabels="c11", values=.2 )
modelACEnone <- omxSetParameters( modelACEnone, labels="c211", newlabels="c11", values=.2 )
modelACEnone <- omxSetParameters( modelACEnone, labels="c311", newlabels="c11", values=.2 )
modelACEnone <- omxSetParameters( modelACEnone, labels="c411", newlabels="c11", values=.2 )
modelACEnone <- omxSetParameters( modelACEnone, labels="c511", newlabels="c11", values=.2 )
modelACEnone <- omxSetParameters( modelACEnone, labels="c611", newlabels="c11", values=.2 )
modelACEnone <- omxSetParameters( modelACEnone, labels="e111", newlabels="e11", values=.3 )
modelACEnone <- omxSetParameters( modelACEnone, labels="e211", newlabels="e11", values=.3 )
modelACEnone <- omxSetParameters( modelACEnone, labels="e311", newlabels="e11", values=.3 )
modelACEnone <- omxSetParameters( modelACEnone, labels="e411", newlabels="e11", values=.3 )
modelACEnone <- omxSetParameters( modelACEnone, labels="e511", newlabels="e11", values=.3 )
modelACEnone <- omxSetParameters( modelACEnone, labels="e611", newlabels="e11", values=.3 )
fitACEnone   <- mxRun( modelACEnone, intervals=T )
sumACEnone   <- summary( fitACEnone )
# mxCompare( fitACErq, fitACE )
mxCompare( fitACEq, fitACEnone )
fitGofs(fitACEnone)
fitEstVCageandsex(fitACEnone)

# Table 2 - columns 2 and 3 (Females and Males)
# Run Sex-Limitation but No Age-Limitation Model
modelACEsex <- mxModel( modelACEq, name="oneACE5sex" )
modelACEsex <- omxSetParameters( modelACEsex, labels="a111", newlabels="a11", values=.6 )
modelACEsex <- omxSetParameters( modelACEsex, labels="a211", newlabels="a11", values=.6 )
modelACEsex <- omxSetParameters( modelACEsex, labels="a311", newlabels="a11", values=.6 )
modelACEsex <- omxSetParameters( modelACEsex, labels="c111", newlabels="c11", values=.2 )
modelACEsex <- omxSetParameters( modelACEsex, labels="c211", newlabels="c11", values=.2 )
modelACEsex <- omxSetParameters( modelACEsex, labels="c311", newlabels="c11", values=.2 )
modelACEsex <- omxSetParameters( modelACEsex, labels="e111", newlabels="e11", values=.3 )
modelACEsex <- omxSetParameters( modelACEsex, labels="e211", newlabels="e11", values=.3 )
modelACEsex <- omxSetParameters( modelACEsex, labels="e311", newlabels="e11", values=.3 )
modelACEsex <- omxSetParameters( modelACEsex, labels="a411", newlabels="a21", values=.6 )
modelACEsex <- omxSetParameters( modelACEsex, labels="a511", newlabels="a21", values=.6 )
modelACEsex <- omxSetParameters( modelACEsex, labels="a611", newlabels="a21", values=.6 )
modelACEsex <- omxSetParameters( modelACEsex, labels="c411", newlabels="c21", values=.2 )
modelACEsex <- omxSetParameters( modelACEsex, labels="c511", newlabels="c21", values=.2 )
modelACEsex <- omxSetParameters( modelACEsex, labels="c611", newlabels="c21", values=.2 )
modelACEsex <- omxSetParameters( modelACEsex, labels="e411", newlabels="e21", values=.3 )
modelACEsex <- omxSetParameters( modelACEsex, labels="e511", newlabels="e21", values=.3 )
modelACEsex <- omxSetParameters( modelACEsex, labels="e611", newlabels="e21", values=.3 )
fitACEsex   <- mxRun( modelACEsex, intervals=T )
sumACEsex   <- summary( fitACEsex )
mxCompare( fitACEq, fitACEsex )
fitGofs(fitACEsex)
fitEstVCageandsex(fitACEsex)

# Table 3
# Run "Adults are different" model with Sex-Limitation Model
modelACEsex2 <- mxModel( modelACEq, name="oneACE5sex2" )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="a111", newlabels="a11", values=.6 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="a211", newlabels="a11", values=.6 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="c111", newlabels="c11", values=.2 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="c211", newlabels="c11", values=.2 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="e111", newlabels="e11", values=.3 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="e211", newlabels="e11", values=.3 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="a411", newlabels="a21", values=.6 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="a511", newlabels="a21", values=.6 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="c411", newlabels="c21", values=.2 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="c511", newlabels="c21", values=.2 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="e411", newlabels="e21", values=.3 )
modelACEsex2 <- omxSetParameters( modelACEsex2, labels="e511", newlabels="e21", values=.3 )
fitACEsex2   <- mxRun( modelACEsex2, intervals=T )
sumACEsex2   <- summary( fitACEsex2 )
# mxCompare( fitACErq, fitACE )
mxCompare( fitACEq, fitACEsex2 )
fitGofs(fitACEsex2)
fitEstVCageandsex(fitACEsex2)